library(Mediana)

# 1. Data model ----
sample_size <- c(100, 150, 200)

# 1.1 Outcome parameter set
# 1.1.1 Standard set
outcome1.control <- parameters(mean = 5.7, sd = 0.30)
outcome1.test <- parameters(mean = 5.6, sd = 0.25)

# 1.1.2 Optimal set
# outcome2.control <- parameter(prop = 0.30)
# outcome2.test <- parameter(prop = 0.50)

# 1.1.3 Pessimeistic set
# outcome3.control <- parameter(dispersion = 0.5, mean = 13)
# outcome3.test <- parameter(dispersion = 0.5, mean = 7.8)

# 1.2 Build Data model
trial.dm <- DataModel() + 
  OutcomeDist(outcome.dist = "NormalDist") +  # Normal distribution endopoints
	SampleSize(sample_size) + 
	Sample(id = "Control",
	       outcome.par = parameters(outcome1.control)) + 
	Sample(id = "Test",
	       outcome.par = parameters(outcome1.test))

# 2. Analysis model ----
# 2.1 Statistic test
trial.am <- AnalysisModel() + 
  Test(id = "Control vs Test",
       samples = samples("Control", "Test"),
       method = "TTest")  # For "Normal distribution endpoint"

# 2.2 Statistic value
trial.am <- trial.am +
	Statistic(id = "Mean_Test",
	          method = "MeanStat",
	          samples = samples("Test"))

# 3. Evaluation model ----
# 3.1 Power analysis
trial.em <- EvaluationModel() + 
  Criterion(id = "Marginal power",
            method = "MarginalPower",
            tests = tests("Control vs Test"),
            labels = c("Control vs Test"),
            par = parameters(alpha = 0.025))

# Statistic value evaluation
trial.em <- trial.em +
	Criterion(id = "Average_Mean",
	          method = "MeanSumm",
	          statistics = statistics("Mean_Test"),
	          labels = c("Average Mean Test"))

# 4. Clinical Scenario Evaluation ----
sim_parameters <- SimParameters(n.sims = 50, proc.load = "full", seed = 4293)
trial.rst <- CSE(trial.dm, trial.am, trial.em, sim_parameters)

# 5. Simulation results ----
# 5.1 Console
trial.rst.summary <- summary(trial.rst)

# 5.2 Presentation model
trial.pm <- PresentationModel() +
	Project(username = "SimTrial",
		title = "Simulation Trial",  # TODO
		description = "Clinical Trial Simulation") +
	Section(by = "outcome.parameter") +
	Table(by = "sample.size") +
	CustomLabel(param = "sample.size",
		label = paste0("N=", sample_size)) +
	CustomLabel(param = "outcome.parameter",
		label = paste0(c("Pessimist", "Standard", "Optimist")))

GenerateReport(
	presentation.model = trial.pm,
	cse.results = trial.rst,
	report.filename = "TrialSimReport.docx")

# 6. Get the data generated in the simulation ----
trial.data <- DataStack(data.model = trial.dm,
                        sim.parameters = sim_parameters)

# i th simulation, j th scenario, k sample
# trial.data$data.set[[i]]$data.scenario[[j]]$sample[[k]]

# Structure of 'sample' (data.frame)
# |-outcome
# |-patient.start.time
# |-patient.end.time
# |-patient.dropout.time
# |-patient.censor.indicator

# Extract simulation results data (data.set)
# sim_data <- ExtractDataStack(data.stack = trial.data,
#                             data.scenario = 1, 
#                             # sample.id = "Control",
#                             simulation.run = 5
#                             )$data.set[[1]]$data.scenario[[1]]

ExtractSimData <- function(trial.data, scenario, run) {
  # Extract simulation results data (data.set)
  sim_data <- ExtractDataStack(data.stack = trial.data,
                               data.scenario = scenario, 
                               # sample.id = "Control",
                               simulation.run = run
  )$data.set[[1]]$data.scenario[[1]]
  
  value1 <- sim_data$sample[[1]]$data$outcome
  value2 <- sim_data$sample[[2]]$data$outcome
  # exchange 0 and 1 (0 -> 1,, 1 -> 0)
  ## 1 means event happened; 0 means no target event meet
  event1 <- 1 - sim_data$sample[[1]]$data$patient.censor.indicator  
  event2 <- 1 - sim_data$sample[[2]]$data$patient.censor.indicator
  id1 <- rep(sim_data$sample[[1]]$id, length(value1))
  id2 <- rep(sim_data$sample[[2]]$id, length(value2))
  
  rst <- data.frame(value = c(value1, value2), id = c(id1, id2),
                    event = c(event1, event2))
  
  return(rst)
}

# 7. Analysis ----

# Simulating Clincail Trial Data
sctd <- ExtractSimData(trial.data, 1, 5)

library(sm)
sm.density.compare(sctd$value, sctd$id, lwd = 2)

boxplot(data ~ id, data = sctd)

library(ggplot2)
# qplot(x = value, data = sctd, geom = "density", color = id, alpha = I(0.5))

ggplot(data = sctd, aes(value, fill = id)) + 
  geom_histogram()
  
ggplot(data = sctd, aes(id, value, color = id)) + 
  geom_boxplot()

ggplot(data = sctd, aes(id, value)) + 
  geom_violin()
